course_project
Bashir,I. (2022)
17/04/2022
Linkage map construction
library(qtl)## Warning: package 'qtl' was built under R version 4.1.3pop3 <- read.cross(format="csv", file="data/population_Z003.csv", genotypes=c("0", "1", "2"), crosstype = "riself")## --Read the following data:
## 191 individuals
## 1106 markers
## 34 phenotypes## Warning in convert2riself(cross): Omitting 12899 genotypes with code==2.## Warning in summary.cross(cross): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
## (Perhaps it is in basepairs?)## --Cross type: riselfmaize <- est.rf(cross = pop3)
plotRF(maize, col.scheme = "redblue")
Marker grouping
max.rf <- 0.38
kosambi <- function(r) (1/4)*log((1+(2*r))/(1-(2*r)))
kosambi(r =max.rf)## [1] 0.4981075For min.lod, we can run Bonferroni correction on the number of tests that we have to perform in order to evaluate marker linkage. The number of tests is the number of marker pairs that we have in our data. As a first guess, we have:
(M <- totmar(maize)) # number of markers## [1] 1106(num.pair <- choose(M,2)) # number of marker pairs## [1] 611065(min.lrt <- qchisq(1-(0.05/num.pair), 1)) # min LRT to consider two markers## [1] 28.7624(min.lod <- 0.2174 * min.lrt) # conversion of LRT to LOD## [1] 6.252946Now, it is time to see how many linkage groups we have. First, we run the function formLinkageGroups() with reorgMarkers = FALSE just to see how the markers get distributed along the formed linkage groups:
maize <- formLinkageGroups(maize, max.rf = 0.38, min.lod = 6.25, reorgMarkers = TRUE)
#table(maize[,2])We note that there are 10 linkage groups (as expected for maize). Now, we can use reorgMarkers = TRUE and update our cross object named maize with each marker numbered according to the linkage group it belongs
plotRF(maize, col.scheme = "redblue")
Marker ordering
We are going to show two ways of ordering markers. The first way uses the orderMarkers() function by R/qtl and usually needs some manual curation. The second way uses the MDS algorithm and it is faster and usually more accurate. You can choose which method you want to use and skip to its specific section, meaning that you do not need to run both ways.
Using orderMarkers function by R/qtl:
R/qtl has a function that runs the Branch-and-Bound algorithm as an possible solution for the Traveling
Salesman Problem (TSP) that ordering markers is. It usually provides a good solution. The problem is that
Branch-and-Bound is very sensitive to the marker choice that is made to start the algorithm. Therefore, we run it at least a couple of times such that some effect of the first choices can be evaluated.
We save the maize object under two new object called maize.bb1 and maize.bb2 , so we can update the
with the results of two runs of the Branch-and-Bound algorithm. In addition, we initialize two objects that will store the log-likelihood of the ordering for each linkage group from both runs:
maize.bb1 <- maize
maize.bb2 <- maize
loglik.bb1 <- loglik.bb2 <- c()Linkage group 1
The object c stores the number of the chromosome under evaluation.
c <- 1
plotRF(maize, chr=c, col.scheme="redblue")
The argument chr let us take a closer look into the heatmap of a specific chromosome, whose markers are clearly not ordered
maize.bb1 <- orderMarkers(cross = maize.bb1, chr=c, use.ripple = FALSE, map.function="kosambi")
plotRF(maize.bb1, chr = c, col.scheme = "redblue")
pull.map(maize.bb1, chr = c)## PHM4913.18 PZA02292.1 PZA03168.5 PZA02737.1 PZA02550.1
## 0.000000e+00 5.000001e-06 5.706188e-01 4.874327e+00 4.874332e+00
## PZA02114.1 PZB01062.3 PZA01476.1 PZA00294.22 PZA03189.4
## 5.449683e+00 6.052846e+00 8.770987e+00 8.770992e+00 8.770997e+00
## PZA01315.1 PZA03561.1 PZA00752.1 PHM5098.25 PZA01267.3
## 9.330274e+00 9.330279e+00 1.935328e+01 2.209257e+01 2.209257e+01
## PZA01135.1 PZA03240.1.2 PZA03465.1 PZA00944.1.2 PZB01235.4
## 2.694937e+01 3.071735e+01 3.168677e+01 3.168677e+01 3.515063e+01
## PHM9418.11 PZA02750.3 PZA02763.1 PZA00939.1 PZA02070.1
## 3.925566e+01 3.925567e+01 3.925567e+01 4.034549e+01 4.241535e+01
## PZA01254.2 PZA02577.1 PZA03200.2 csu1138.3.4 PZA02741.1
## 4.241535e+01 4.308020e+01 4.308020e+01 4.308021e+01 5.183141e+01
## an1.5 PZA02135.2 PZA00455.14.16 PHM1968.22 PZA02191.1
## 5.183141e+01 5.183142e+01 5.349319e+01 5.532299e+01 5.532300e+01
## PZA00068.1 PZA03531.1 PZA00619.3 PZA02467.10 PZA03074.27
## 5.632409e+01 5.662325e+01 7.055372e+01 7.146288e+01 7.856654e+01
## PZA01391.1 PZA01216.1 PZA03194.1 PZA01019.1 PZA01963.15
## 7.948514e+01 7.948515e+01 7.948515e+01 7.948516e+01 7.948516e+01
## PZA00131.15 PHM5480.17 PZA03193.2 PHM12706.14 PZA01039.1
## 7.948517e+01 7.948517e+01 8.649254e+01 8.716110e+01 8.716519e+01
## PZA02014.3 PZA03265.3 PZA03741.1 PHM6043.19 PHM15871.11
## 8.716519e+01 8.908646e+01 9.105424e+01 9.105424e+01 9.410745e+01
## PHM5484.22 PZA02117.1 PZA00658.21 PHM2478.22 PZA02823.1
## 9.410746e+01 9.410746e+01 9.690109e+01 9.751182e+01 9.751183e+01
## umc128.2 PHM4942.12 PZA00664.3 PZA02186.1 PZB01647.1
## 9.843189e+01 9.843189e+01 9.843190e+01 9.873443e+01 9.993908e+01
## PZA03001.15 PZA00381.4 PZA03301.2 PHM4926.16 PZA03064.6
## 9.993908e+01 1.015890e+02 1.030467e+02 1.030467e+02 1.068317e+02
## PHM16605.19 PZA02269.3.4 PZA03404.1 kip1.3 PHM3034.3
## 1.081519e+02 1.132823e+02 1.132824e+02 1.154321e+02 1.154321e+02
## PHM14475.7 glb1.2 PZA01588.1 PZA01921.20.19 PHM5526.25
## 1.189909e+02 1.189909e+02 1.189909e+02 1.195877e+02 1.195877e+02
## PZA00339.4 PZA02985.5 PZA03457.1 PZB00063.1 PZB00895.1
## 1.195877e+02 1.210303e+02 1.213131e+02 1.261977e+02 1.261977e+02
## PZB00008.1 PZA02278.1 PZA02698.3 PZA02520.1 PZA00030.11
## 1.261978e+02 1.261978e+02 1.283516e+02 1.296424e+02 1.296424e+02
## PZB00114.1 PZA03188.3 PZA02204.1 PZA01978.23 PZA00610.16
## 1.306897e+02 1.378545e+02 1.378545e+02 1.378545e+02 1.378545e+02
## PZA03020.8 PZA02957.5 PZA00978.1 PZA02087.2 PZA01246.1
## 1.382621e+02 1.382621e+02 1.394943e+02 1.408175e+02 1.408175e+02
## PZA00245.20 PHM18705.23 PZB01403.1 PZA03037.2 PZA03305.7.1
## 1.419242e+02 1.422860e+02 1.426478e+02 1.441362e+02 1.441362e+02
## PZB01227.6 PZA00894.7 PZA02044.1 PZA00307.14 PZA00991.2
## 1.441362e+02 1.441362e+02 1.510183e+02 1.549186e+02 1.569589e+02
## PZA00235.9 PZA00276.18 PZA00623.3 PZA02359.10 PHM9807.9
## 1.569589e+02 1.569589e+02 1.585866e+02 1.585866e+02 1.609158e+02
## PZA01238.1.2 PZA00343.31 PZA01068.1 PHM1275.22 PZA00243.25
## 1.621469e+02 1.624272e+02 1.624272e+02 1.627476e+02 1.641581e+02
## PZA00856.2 PHM7616.35 PZA01239.2 PZA01807.1 PZA00432.4
## 1.641581e+02 1.681775e+02 1.681775e+02 1.681775e+02 1.681775e+02
## PZA03613.1 PZA01271.1 PZA02129.1 PZA02032.1 PHM2244.142
## 2.411258e+02 2.414039e+02 2.434352e+02 2.440836e+02 2.456885e+02
## PHM6238.36 PZA00181.2 PZA02372.1 PZA00528.1 PZA00447.8
## 2.504676e+02 2.504676e+02 2.504676e+02 2.533944e+02 2.537385e+02
## PZA00175.2 PZA00566.5 PZA02284.1 PZA00731.7 PZA03521.1
## 2.537385e+02 2.558706e+02 2.558706e+02 2.558706e+02 2.582544e+02
## PZA00106.10 PZA00887.1 csu1171.2 PZA03551.1 PZA01497.1
## 2.582544e+02 2.624199e+02 2.705287e+02 2.705287e+02 2.705287e+02
## PZA02094.9 PZA01652.1 PZA02393.2 PZB00718.5 PZB00648.5
## 2.705287e+02 2.705287e+02 2.745808e+02 2.759369e+02 2.759369e+02
## PZA01030.1 PHM3226.15 PZA02487.1 PZA00425.11 PHM13619.5
## 2.762315e+02 2.844797e+02 2.844797e+02 2.844797e+02 2.844797e+02
## PHM4531.46 PZB01957.1 PZB02058.1 PZA02686.1 PZA01455.1
## 2.906317e+02 2.906317e+02 2.935355e+02 3.032138e+02 3.032138e+02
## PZA02490.1 PZB01662.1 PZA01348.1 PZA02271.1 PZA02195.1
## 3.032138e+02 3.032138e+02 3.032138e+02 3.046707e+02 3.046707e+02
## PZA02376.1 PHM3726.129 PZA00240.6 PZA00081.18 PZA00962.1
## 3.080497e+02 3.080497e+02 3.080497e+02 3.089379e+02 3.089379e+02
## PZA03243.2 PZA03742.1 PZA03183.5 umc13.1 PZB00872.3
## 3.089379e+02 3.089379e+02 3.098421e+02 3.098421e+02 3.104436e+02(loglik.bb1[c] <- attr(maize.bb1$geno[[c]]$map,"loglik"))## [1] -3032.742maize.bb2 <- orderMarkers(cross = maize.bb2, chr=c, use.ripple = FALSE, map.function="kosambi")
plotRF(maize.bb2, chr = c, col.scheme = "redblue")
pull.map(maize.bb2, chr = c)## PZB00872.3 PZA03183.5 umc13.1 PZA00962.1 PZA03243.2
## 0.0000000 0.6014728 0.6014778 1.5057287 1.5057337
## PZA03742.1 PZA02376.1 PZA00081.18 PZA02195.1 PHM3726.129
## 1.5057387 1.5057437 1.5057487 2.3939367 2.3939417
## PZA00240.6 PZA02271.1 PZA01455.1 PZA02686.1 PZA01348.1
## 2.3939467 5.7732633 7.2298510 7.2298560 7.2298610
## PZB01662.1 PZB02058.1 PZA02490.1 PZB01957.1 PZA00425.11
## 7.2298660 16.9072498 16.9072548 19.8109630 25.9637887
## PZA02487.1 PHM13619.5 PHM4531.46 PHM3226.15 PZA01030.1
## 25.9637937 25.9637987 25.9638037 25.9638087 34.2119525
## PZB00648.5 PZB00718.5 PZA02393.2 PZA01497.1 PZA03551.1
## 34.5065893 34.5065943 35.8627121 39.9147296 39.9147346
## PZA01652.1 csu1171.2 PZA02094.9 PZA00887.1 PZA00566.5
## 39.9147396 39.9147446 39.9147496 48.0235816 52.1890624
## PZA03521.1 PZA00106.10 PZA02284.1 PZA00731.7 PZA00175.2
## 52.1890674 52.1890724 54.5729129 54.5729179 56.7050100
## PZA00447.8 PZA00528.1 PZA00181.2 PHM6238.36 PHM2244.142
## 56.7050150 57.0490959 59.9755998 59.9756048 64.7545038
## PZA02372.1 PZA02032.1 PZA02129.1 PZA01271.1 PZA03613.1
## 64.7545088 66.3597717 67.0081745 69.0395631 69.3176620
## PZA00432.4 PHM7616.35 PZA01807.1 PZA01239.2 PZA00856.2
## 142.1778828 142.1778878 142.1778928 142.1778978 146.1971823
## PZA00243.25 PHM1275.22 PZA01068.1 PZA00343.31 PZA01238.1.2
## 146.1971873 147.6076243 147.9280489 147.9280539 148.2088931
## PHM9807.9 PZA02359.10 PZA00623.3 PZA00991.2 PZA00307.14
## 149.4398011 151.7774712 151.7774762 153.4109356 155.0707978
## PZA00235.9 PZA00276.18 PZA02044.1 PZA00894.7 PZA03305.7.1
## 156.5436236 156.5436286 162.2192945 168.7978964 168.7979014
## PZB01227.6 PZA03037.2 PZB01403.1 PHM18705.23 PZA00245.20
## 168.7979064 168.7979114 170.2510736 170.6147530 170.9784657
## PZA02087.2 PZA00978.1 PZA01246.1 PZA03188.3 PZA03020.8
## 172.0910165 173.4218023 173.4218073 174.6568515 174.6568565
## PZA02957.5 PZA02204.1 PZA01978.23 PZA00610.16 PZB00114.1
## 174.6568615 175.0649641 175.0649691 175.0649741 182.2356230
## PZA02520.1 PZA00030.11 PZA02698.3 PZB00008.1 PZB00063.1
## 183.2831073 183.2831123 184.5741989 186.7281796 186.7281846
## PZA02278.1 PZB00895.1 PZA03457.1 PZA02985.5 PZA01921.20.19
## 186.7281896 186.7281946 191.6133768 191.8961806 193.3389500
## PHM5526.25 PZA00339.4 glb1.2 PZA01588.1 PHM3034.3
## 193.3389550 193.3435091 193.9356389 193.9356439 197.4944569
## kip1.3 PHM14475.7 PZA02269.3.4 PHM16605.19 PZA03404.1
## 197.4944619 197.4944669 199.6441497 204.7745574 204.7745624
## PZA03064.6 PHM4926.16 PZA03301.2 PZA00381.4 PZA03001.15
## 206.0952684 209.8801478 209.8801528 211.3378762 212.9878000
## PZB01647.1 PZA02186.1 PHM4942.12 umc128.2 PZA02823.1
## 212.9878050 214.1924517 214.4949862 214.4949912 215.4150041
## PHM2478.22 PZA00664.3 PZA00658.21 PZA03741.1 PHM15871.11
## 215.4150091 215.4150141 216.0257951 218.8194143 218.8194193
## PZA02117.1 PHM5484.22 PHM6043.19 PZA03265.3 PZA02014.3
## 218.8194243 218.8194293 221.8729957 223.8404649 225.7617373
## PZA01039.1 PZA03193.2 PZA03194.1 PHM12706.14 PZA01391.1
## 225.7617423 226.4351550 226.4351600 226.4351650 233.4416391
## PZA00131.15 PHM5480.17 PZA01216.1 PZA01963.15 PZA01019.1
## 233.4416441 233.4416491 233.4416541 233.4416591 233.4416641
## PZA03074.27 PZA02467.10 PZA00619.3 PZA03531.1 PHM1968.22
## 234.3602674 241.4639129 242.3730707 256.3034845 256.6027474
## PZA00068.1 PZA02191.1 PZA00455.14.16 an1.5 csu1138.3.4
## 256.6027524 257.6038290 259.4335309 261.0944026 269.8469246
## PZA02741.1 PZA03200.2 PZA02577.1 PZA02135.2 PZA01254.2
## 269.8469296 269.8469346 269.8469396 269.8469446 270.5117869
## PZA02070.1 PZA00939.1 PZA02750.3 PZA02763.1 PHM9418.11
## 270.5117919 272.5817120 273.6715603 273.6715653 273.6715703
## PZB01235.4 PZA00944.1.2 PZA03465.1 PZA03240.1.2 PZA01135.1
## 277.7767056 281.2414377 281.2414427 282.2110969 285.9824701
## PZA01267.3 PZA00752.1 PZA01315.1 PZA03561.1 PHM5098.25
## 290.7836408 293.4755324 303.4990947 303.4990997 304.0583705
## PZA00294.22 PZA03189.4 PZA01476.1 PZB01062.3 PZA02114.1
## 304.0583755 304.0583805 304.0583855 306.7765075 307.3794684
## PZA02550.1 PZA02292.1 PHM4913.18 PZA02737.1 PZA03168.5
## 307.9545896 312.2920840 312.2920890 312.8787447 312.8787497Linkage group 2
c <- 2
plotRF(maize, chr=c, col.scheme = "redblue")
maize.bb1 <- orderMarkers(cross = maize.bb1, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb1, chr = c)## PZA01060.1 PZA02769.1 PZA02480.1 PZA02390.1 PZA01140.1
## 0.000000e+00 5.000001e-06 1.521297e+00 1.611027e+01 1.675923e+01
## PZA00836.1 PZA01259.1 PZA01680.3 PZA03167.5 PZA02015.11
## 1.675948e+01 2.464741e+01 2.469483e+01 2.588630e+01 2.735286e+01
## PZA00545.26 PZA02068.1 PZA00980.1 PZA00963.3 PHM3512.186
## 2.935686e+01 3.601359e+01 3.751027e+01 3.784517e+01 3.959568e+01
## PZA00395.2 PZA02667.1 PZB00765.1 PZA01265.1 PZA02513.1
## 4.037148e+01 4.154840e+01 4.219372e+01 4.219372e+01 4.219373e+01
## PZA02060.1 PZA02820.17 PZA01142.4 PZA00652.17 PZA03024.16
## 4.219373e+01 4.362748e+01 4.506669e+01 4.539625e+01 4.539625e+01
## PHM532.23 PZA02411.3 PZA02426.1 PZA00352.23 PZA02751.1
## 4.996339e+01 5.076180e+01 5.156058e+01 5.250475e+01 5.343755e+01
## PZA03320.6 PZA03317.1 PZA03172.3 PZA02383.1 PZB01017.1
## 5.343755e+01 5.405840e+01 5.435953e+01 5.435953e+01 7.604198e+01
## PZA01608.1 PZA01796.1 ae1.8.7 PZA01763.2 PZA02981.2
## 7.604199e+01 7.731341e+01 8.035158e+01 8.035159e+01 8.035159e+01
## PZA00148.3 PZA01294.2.1 PZA03536.1 PZA02641.2 PZA00255.14
## 8.035160e+01 8.146505e+01 8.146505e+01 8.146506e+01 8.177995e+01
## PZA00300.14 PZA01410.1 PZA00987.1 PZA02040.2 PZA03717.1
## 8.258298e+01 8.312694e+01 8.312694e+01 8.358165e+01 8.641484e+01
## PHM1899.157 PZA01304.1 PZA03714.1 PZA02633.4 PZA02209.2
## 8.786312e+01 8.786312e+01 8.786313e+01 8.852906e+01 8.889514e+01
## PZA02408.2 PZA02356.7 PHM5296.6 PZA03324.1 PZA01575.1
## 8.889514e+01 8.889515e+01 8.889515e+01 8.889516e+01 8.962350e+01
## PZA03452.6 PZA00067.10 PZA01365.1 PZA02164.16 PZA00881.1
## 9.113211e+01 1.105422e+02 1.126569e+02 1.137952e+02 1.149271e+02
## PZA00643.13 PZA03049.24 PZA01693.1 PZA00273.5 PZA01779.1
## 1.154709e+02 1.157661e+02 1.170803e+02 1.170803e+02 1.183945e+02
## PZA02818.6 PZA02862.3 PZA00261.6 PHM5798.39 PZA02525.1
## 1.188964e+02 1.194682e+02 1.194682e+02 1.209574e+02 1.209574e+02
## PZA01303.1 PZA03677.1 PZA01349.2 PZB01112.1 PHM4165.14
## 1.209574e+02 1.209574e+02 1.213216e+02 1.216856e+02 1.227988e+02
## PZA01050.1 PZB00232.2 PZA03451.5 PZA02676.2 PHM1870.20
## 1.231626e+02 1.231626e+02 1.234753e+02 1.234753e+02 1.234753e+02
## PHM3171.5 PZB01115.3 PZA00222.7 PZA00522.12.7 PZA01804.1
## 1.237621e+02 1.237621e+02 1.252321e+02 1.255188e+02 1.255188e+02
## PZA00805.1 PZA02207.1 PZA00981.3 PHM12992.5 PZA00996.1
## 1.258429e+02 1.288016e+02 1.315967e+02 1.315968e+02 1.315968e+02
## PHM3691.18 PZA01530.1 PHM16854.3 PHM4647.8 PZA00499.3
## 1.315968e+02 1.315968e+02 1.319554e+02 1.319554e+02 1.319554e+02
## PZA01563.1 PZA00801.1 PZB00869.4 PZA00934.2 PZA02113.1
## 1.319554e+02 1.323241e+02 1.323241e+02 1.345911e+02 1.345912e+02
## PHM565.31 PZA01427.1 PZA03226.3 PZA02792.26.25 PZA03274.4
## 1.351010e+02 1.351010e+02 1.401091e+02 1.403524e+02 1.403524e+02
## PZA00517.7 PZA01523.1 PZA03578.1 PZA01327.1 PZA00985.1
## 1.422119e+02 1.436708e+02 1.441564e+02 1.451431e+02 1.484146e+02
## PZA00112.5 PZA03092.7 PZA01284.6 PZB00079.4 PZA01371.1
## 1.501857e+02 1.523747e+02 1.545730e+02 1.567795e+02 1.610688e+02
## PZA01925.1 PZA00865.1 PZA02029.21 PHM3137.17 PZB00094.1
## 1.610688e+02 1.610688e+02 1.610688e+02 1.613730e+02 1.613730e+02
## PZA02753.1 PZA02462.1 PZB00054.3 PZA02653.12 PHM13122.43
## 1.634001e+02 1.646080e+02 1.646080e+02 1.694607e+02 1.766544e+02
## PHM5359.10 PZA01570.1 PZA02316.22 PZA01438.1 PZA00191.5
## 1.778748e+02 1.778748e+02 1.800035e+02 1.813322e+02 1.848231e+02
## PZA00818.1 PZA01983.1 PZA02367.1 PZA01887.1
## 1.897547e+02 1.897547e+02 1.897547e+02 1.897547e+02(loglik.bb1[c] <- attr(maize.bb1$geno[[c]]$map, "loglik"))## [1] -2304.4maize.bb2 <- orderMarkers(cross = maize.bb2, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb2, chr = c)## PZA01060.1 PZA02769.1 PZA02480.1 PZA02390.1 PZA01140.1
## 0.000000e+00 5.000001e-06 1.521173e+00 1.612776e+01 1.677974e+01
## PZA01680.3 PZA00836.1 PZA01259.1 PZA03167.5 PZA02015.11
## 2.431275e+01 2.431276e+01 2.463947e+01 2.571958e+01 2.717412e+01
## PZA00545.26 PZA02068.1 PZA00980.1 PZA00963.3 PHM3512.186
## 2.918244e+01 3.583325e+01 3.732967e+01 3.766450e+01 3.941462e+01
## PZA00395.2 PZA02667.1 PZA01265.1 PZA02060.1 PZB00765.1
## 4.019043e+01 4.136736e+01 4.201281e+01 4.201293e+01 4.201293e+01
## PZA02513.1 PZA02820.17 PZA01142.4 PZA03024.16 PZA00652.17
## 4.201294e+01 4.344684e+01 4.488612e+01 4.521572e+01 4.521573e+01
## PHM532.23 PZA02411.3 PZA02426.1 PZA03320.6 PZA00352.23
## 4.978066e+01 5.057878e+01 5.137717e+01 5.329389e+01 5.329390e+01
## PZA02751.1 PZA03172.3 PZA02383.1 PZA03317.1 PZB01017.1
## 5.329390e+01 5.427153e+01 5.459979e+01 5.459980e+01 7.560648e+01
## PZA01608.1 PZA01796.1 PZA01763.2 ae1.8.7 PZA02981.2
## 7.560648e+01 7.687853e+01 7.992167e+01 7.992167e+01 7.992168e+01
## PZA00148.3 PZA01294.2.1 PZA03536.1 PZA02641.2 PZA00255.14
## 7.992168e+01 8.103699e+01 8.103699e+01 8.103700e+01 8.135232e+01
## PZA02040.2 PZA00300.14 PZA01410.1 PZA00987.1 PZA03717.1
## 8.215611e+01 8.215612e+01 8.270064e+01 8.270065e+01 8.499160e+01
## PHM1899.157 PZA03714.1 PZA02209.2 PZA03324.1 PZA02356.7
## 8.645594e+01 8.645595e+01 8.758093e+01 8.758094e+01 8.758094e+01
## PZA03452.6 PZA01575.1 PZA02408.2 PHM5296.6 PZA02633.4
## 8.962543e+01 9.097170e+01 9.169378e+01 9.169379e+01 9.200661e+01
## PZA01304.1 PZA00067.10 PZA02164.16 PZA01365.1 PZA00881.1
## 9.257205e+01 1.084400e+02 1.105809e+02 1.105809e+02 1.128856e+02
## PZA00643.13 PZA03049.24 PZA01693.1 PZA00273.5 PZA01779.1
## 1.134223e+02 1.137137e+02 1.150100e+02 1.150100e+02 1.163063e+02
## PZA02818.6 PZA00261.6 PZA02862.3 PZA01303.1 PZA03677.1
## 1.168027e+02 1.173675e+02 1.173675e+02 1.188461e+02 1.188461e+02
## PZA02525.1 PHM5798.39 PZA01349.2 PZB01112.1 PHM4165.14
## 1.188461e+02 1.188461e+02 1.192117e+02 1.195772e+02 1.206952e+02
## PHM3171.5 PZA02676.2 PZA03451.5 PZB01115.3 PZB00232.2
## 1.210604e+02 1.210604e+02 1.210604e+02 1.210604e+02 1.210604e+02
## PZA01050.1 PHM1870.20 PZA00222.7 PZA01804.1 PZA00522.12.7
## 1.210604e+02 1.213483e+02 1.225219e+02 1.228098e+02 1.228098e+02
## PZA00805.1 PZA02207.1 PZA01530.1 PHM3691.18 PZA00996.1
## 1.231354e+02 1.260973e+02 1.288958e+02 1.288958e+02 1.288959e+02
## PZA00981.3 PHM12992.5 PHM16854.3 PHM4647.8 PZA01563.1
## 1.288959e+02 1.288959e+02 1.292552e+02 1.292552e+02 1.292552e+02
## PZA00499.3 PZA00801.1 PZB00869.4 PZA02113.1 PZA00934.2
## 1.292552e+02 1.296258e+02 1.296258e+02 1.307306e+02 1.318378e+02
## PHM565.31 PZA02792.26.25 PZA01427.1 PZA03226.3 PZA03274.4
## 1.323542e+02 1.323542e+02 1.323542e+02 1.373857e+02 1.374815e+02
## PZA00517.7 PZA01523.1 PZA03578.1 PZA01327.1 PZA00985.1
## 1.391848e+02 1.406449e+02 1.411304e+02 1.421171e+02 1.453851e+02
## PZA03092.7 PZA00112.5 PZA01284.6 PZB00079.4 PZA02029.21
## 1.471551e+02 1.471551e+02 1.517683e+02 1.539661e+02 1.582528e+02
## PZA01371.1 PZA01925.1 PZA00865.1 PHM3137.17 PZB00094.1
## 1.582528e+02 1.582528e+02 1.582528e+02 1.585569e+02 1.585569e+02
## PZA02753.1 PZB00054.3 PZA02653.12 PZA02462.1 PHM13122.43
## 1.605833e+02 1.617909e+02 1.666435e+02 1.666435e+02 1.738392e+02
## PZA01570.1 PHM5359.10 PZA02316.22 PZA01438.1 PZA00191.5
## 1.750600e+02 1.761004e+02 1.771406e+02 1.784699e+02 1.819620e+02
## PZA01983.1 PZA01887.1 PZA02367.1 PZA00818.1
## 1.868940e+02 1.868940e+02 1.868940e+02 1.868940e+02plotRF(maize.bb2, chr = c, col.scheme = "redblue")
(loglik.bb2[c] <- attr(maize.bb2$geno[[c]]$map, "loglik"))## [1] -2324.35Linkage group 3
c <- 3
plotRF(maize, chr=c, col.scheme = "redblue")
maize.bb1 <- orderMarkers(cross = maize.bb1, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb1, chr = c)## PZA00309.1 PZA02090.1 PZD00038.2 PZA02678.1 PHM12859.7
## 0.00000 13.46322 16.47789 18.06638 18.36297
## PZA00100.10 PZA03527.1 PZA03212.3 PZA00749.1 PZB01944.1
## 18.36297 20.92875 25.17167 27.48102 29.79103
## PZA01765.1 PZA00508.2 PZA02098.2 PHM4204.69 PZA02427.1
## 30.69109 30.69110 30.69110 38.80000 47.71134
## PHM2343.25 PHM4145.18 PZA03054.5 PZA01473.1 PZA02255.2
## 47.71134 47.71135 54.28889 54.28889 54.28890
## PZA00348.11 PZA00210.1.9 PHM13823.7 PZA01114.2 PZA00297.2
## 54.28890 54.28891 54.96652 54.96652 54.96653
## zb21.1 PZA00380.10 PHM15899.9 PZA03070.9 PZA00265.6
## 54.96653 54.96654 55.42202 55.42203 57.40815
## PZA00279.2 PHM15474.5 PZA03119.1 PZA00581.3 PZA02589.1
## 57.40816 57.40818 57.40818 57.40819 57.40819
## PZA01447.1 PZA02296.1 PZA02699.1 PZA00509.1 PZA02134.3
## 57.40820 58.92027 58.92028 58.92028 59.27384
## PHM5502.31 PZA02742.1 PZA00707.9 PZA02619.1 PZA02645.2
## 59.62733 59.98732 60.80094 61.61320 61.61320
## PZA03198.3 PHM15449.10 PZA00363.7 PZA00413.20.18 PZD00016.4
## 63.47126 64.32849 64.32849 64.32850 66.82813
## PZA02299.16 PZB02002.1 PZD00015.5 PZA02474.1 PHM1745.16
## 66.82814 66.82814 66.82815 68.79768 68.79769
## PZA00920.1 PHM890.20 PZB02044.1 PHM4955.12 PZB02122.1
## 69.13376 70.16722 70.44483 71.79349 71.79350
## PZA01934.6 PZA00827.1 PZA00948.1 PHM9914.11 PZB02179.1
## 71.79350 72.98389 74.19445 74.47859 74.47859
## PZA00828.2 PZA02423.1 PZA00088.3 PHM2423.33 PZA02182.1
## 74.47860 124.91680 125.47299 128.29801 134.81914
## PZA01360.3 PZA01688.3 PZA00316.10 PZA00402.1 PHM3852.23
## 135.67582 136.54456 136.54457 136.54457 136.54458
## PZA00219.7 PHM2672.19 PZA03391.1 PZA02668.2 PZA01154.1
## 148.37946 148.37947 148.37947 149.44296 151.40824
## PHM3342.31 PZA01233.1 PZA02824.4 PZA02514.1 PZA00750.1
## 151.40824 151.40825 151.40825 151.40826 154.47061
## PZA02665.2 sh2.21 PZA03146.4 PZB01457.1 PZA02616.1
## 154.47062 156.55749 157.22980 161.19239 164.92586
## PZA02516.1 PZA00538.18.15 PZA02733.1 PZA03154.4 PZA00892.5
## 167.26084 168.84410 173.63479 176.06518 177.32858
## PZA01501.1 PZA02122.9 PZA01035.1 PZB01109.1 PZA01457.1
## 178.97303 179.99265 182.66121 182.66122 183.03695
## PZA00308.24 PZA03255.1 PHM13673.53 PZA03743.1 PZA03744.1
## 183.03702 184.88202 185.19368 185.50550 185.50550
## PZA01228.2 PHM824.17 PZA00494.2 PZA03647.1 PZA03191.1.4
## 185.81731 185.81732 187.48954 188.78123 189.80811
## zb27.1 PZA03735.1 PZA03733.1 PHM1675.29 PZA02212.1
## 189.80811 194.01180 194.01180 194.29656 197.60866
## PZA02654.3 PHM17210.5 PZA01962.12 PZA01726.1 PZA00783.1
## 197.60866 200.12257 202.63679 205.15488 206.82343
## PZA02402.1 PZA03032.19 PZA03073.28.26 PHM1959.26 PZD00027.2
## 208.22774 208.22774 209.23787 209.23787 209.23788
## PHM2885.31 PZA01396.1 PHM4621.57 PZA00186.4 PZA00667.2
## 210.25333 210.95397 210.95398 210.95398 214.36524plotRF(maize.bb1, chr = c, col.scheme = "redblue")
(loglik.bb1[c] <- attr(maize.bb1$geno[[c]]$map, "loglik"))## [1] -2249.112maize.bb2 <- orderMarkers(cross = maize.bb2, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb2, chr = c)## PZA00309.1 PZA02090.1 PZD00038.2 PZA02678.1 PZA00100.10
## 0.00000 13.46198 16.47562 18.06270 18.35889
## PZA03527.1 PHM12859.7 PZA00749.1 PZA03212.3 PZB01944.1
## 20.92004 20.92005 25.12571 25.12571 30.02652
## PZA02098.2 PHM4204.69 PZA01765.1 PZA00508.2 PZA00210.1.9
## 30.92504 39.02418 39.02419 39.02419 47.93946
## PHM2343.25 PHM4145.18 PZA03054.5 PZA02427.1 PZA02255.2
## 47.93947 47.93947 54.50797 54.50798 54.50798
## PZA01473.1 PZA00348.11 PHM13823.7 PZA01114.2 PZA00380.10
## 54.50799 54.50799 55.18897 55.18898 55.18898
## zb21.1 PZA00297.2 PHM15899.9 PZA00581.3 PZA03070.9
## 55.18899 55.18899 55.63027 55.63027 55.63028
## PZA02589.1 PZA00279.2 PZA00265.6 PHM15474.5 PZA01447.1
## 57.46495 57.46499 57.46502 57.46506 57.46506
## PZA03119.1 PZA00509.1 PZA02699.1 PZA02296.1 PZA02134.3
## 57.77699 59.53288 59.53289 59.53289 59.88272
## PHM5502.31 PZA00707.9 PZA02742.1 PZA02645.2 PZA02619.1
## 60.23247 60.59087 60.59087 61.39911 62.20869
## PZA03198.3 PZA00413.20.18 PHM15449.10 PZD00015.5 PZA02299.16
## 64.06313 64.92041 64.92042 67.42025 67.42026
## PZB02002.1 PZA00363.7 PZD00016.4 PZA02474.1 PHM1745.16
## 67.42026 67.42027 67.42027 69.38986 69.38986
## PZA00920.1 PHM890.20 PZB02044.1 PHM4955.12 PZA01934.6
## 69.72594 70.75940 71.03701 72.38568 72.38569
## PZB02122.1 PZA00827.1 PZA00948.1 PHM9914.11 PZB02179.1
## 72.38569 73.57609 74.78667 75.07081 75.07081
## PZA00828.2 PZA02423.1 PZA00088.3 PHM2423.33 PZA02182.1
## 75.07082 125.51493 126.07112 128.89615 135.41729
## PZA01360.3 PZA00316.10 PZA01688.3 PHM3852.23 PHM2672.19
## 136.27397 137.14273 137.14274 137.14274 148.97684
## PZA00402.1 PZA03391.1 PZA00219.7 PZA02668.2 PHM3342.31
## 148.97685 148.97685 148.97686 150.04165 152.00851
## PZA02665.2 PZA01233.1 PZA02514.1 PZA02824.4 PZA01154.1
## 152.00852 152.00852 152.00853 152.00853 153.48786
## PZA00750.1 sh2.21 PZA03146.4 PZB01457.1 PZA02616.1
## 154.97003 157.06118 157.73482 161.69822 165.43273
## PZA02516.1 PZA00538.18.15 PZA02733.1 PZA03154.4 PZA00892.5
## 167.76799 169.35122 174.14129 176.57294 177.84490
## PZA01501.1 PZA02122.9 PZA01457.1 PZA01035.1 PZB01109.1
## 179.50411 180.56927 182.42601 183.15406 183.15406
## PZA00308.24 PZA03255.1 PHM13673.53 PZA03743.1 PZA03744.1
## 183.36239 185.10315 185.41491 185.72681 185.72681
## PHM824.17 PZA01228.2 PZA00494.2 PZA03647.1 zb27.1
## 186.03871 186.03872 187.71032 189.00225 190.02932
## PZA03191.1.4 PZA03733.1 PZA03735.1 PHM1675.29 PZA01962.12
## 190.02932 194.23397 194.23398 194.51872 197.83056
## PZA02654.3 PZA02212.1 PHM17210.5 PZA01726.1 PZA03032.19
## 197.83057 200.34462 202.85891 205.37625 207.04492
## PZA00783.1 PZA02402.1 PZA03073.28.26 PHM1959.26 PZD00027.2
## 207.04492 208.44923 209.45922 209.45923 209.45923
## PHM2885.31 PHM4621.57 PZA01396.1 PZA00186.4 PZA00667.2
## 210.47468 211.17533 211.17533 211.17534 214.58660(loglik.bb2[c] <- attr(maize.bb2$geno[[c]]$map, "loglik"))## [1] -2243.235Linkage group 4
for(i in 1:10) {
c <- i
plotRF(maize, chr=c, col.scheme = "redblue")
maize.bb1 <- orderMarkers(cross = maize.bb1, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb1, chr = c)
plotRF(maize.bb1, chr = c, col.scheme = "redblue")
(loglik.bb1[c] <- attr(maize.bb1$geno[[c]]$map, "loglik"))
maize.bb2 <- orderMarkers(cross = maize.bb2, chr=c, use.ripple = FALSE, map.function="kosambi")
pull.map(maize.bb2, chr = c)
(loglik.bb2[c] <- attr(maize.bb2$geno[[c]]$map, "loglik"))
}
Manual adjustments from orderMarkers() output
Before any manual adjustments to the marker orderings, we need to check which runs of orderMarkers()performed best for each linkage group. In order to do so, we should look at linkage group lengths, maximum space between two consecutive markers (a.k.a. gap) and their log-likelihoods. Higher log-likelihoods indicate better (i.e. more likely) maps.
We combine the information from the function summaryMap() with the one gathered in the objects loglik.bb1 and loglik.bb1 from each respective run:
knitr::kable(cbind(summaryMap(maize.bb1), log.likelihood=c(loglik.bb1, sum(loglik.bb1))))| n.mar | length | ave.spacing | max.spacing | log.likelihood | |
|---|---|---|---|---|---|
| 1 | 175 | 281.2302 | 1.616265 | 46.61441 | -3047.048 |
| 2 | 139 | 176.7643 | 1.280901 | 18.12182 | -2312.367 |
| 3 | 130 | 229.0488 | 1.775572 | 57.14253 | -2309.444 |
| 4 | 127 | 281.1182 | 2.231096 | 84.08295 | -2347.497 |
| 5 | 111 | 203.5704 | 1.850640 | 54.51587 | -2127.110 |
| 6 | 106 | 184.8542 | 1.760516 | 56.07134 | -1750.309 |
| 7 | 85 | 194.2696 | 2.312733 | 84.18438 | -1555.094 |
| 8 | 78 | 195.3088 | 2.536478 | 78.35844 | -1522.708 |
| 9 | 78 | 265.3874 | 3.446589 | 88.10087 | -1679.812 |
| 10 | 77 | 196.8446 | 2.590061 | 79.54214 | -1547.089 |
| overall | 1106 | 2208.3963 | 2.014960 | 88.10087 | -20198.478 |
plotMap(maize.bb1)
plotRF(maize.bb1, col.scheme = "redblue")
knitr::kable(cbind(summaryMap(maize.bb2), log.likelihood=c(loglik.bb2, sum(loglik.bb2))))| n.mar | length | ave.spacing | max.spacing | log.likelihood | |
|---|---|---|---|---|---|
| 1 | 175 | 305.0167 | 1.752970 | 72.94825 | -2992.181 |
| 2 | 139 | 176.4208 | 1.278412 | 18.03897 | -2289.328 |
| 3 | 130 | 232.4914 | 1.802259 | 57.16268 | -2344.634 |
| 4 | 127 | 274.1921 | 2.176128 | 69.34898 | -2336.042 |
| 5 | 111 | 203.5703 | 1.850639 | 54.51589 | -2127.111 |
| 6 | 106 | 185.1711 | 1.763534 | 56.06773 | -1749.789 |
| 7 | 85 | 213.9895 | 2.547494 | 84.20014 | -1584.014 |
| 8 | 78 | 117.4314 | 1.525083 | 10.87745 | -1406.901 |
| 9 | 78 | 211.5313 | 2.747159 | 88.08805 | -1545.376 |
| 10 | 77 | 166.5984 | 2.192084 | 61.31265 | -1474.740 |
| overall | 1106 | 2086.4130 | 1.903661 | 88.08805 | -19850.117 |
plotMap(maize.bb2)
plotRF(maize.bb2, col.scheme = "redblue")
save.image("maize_z003.RData")Now, we can select the best order so far, that we can still try to improve by making manual adjustments for each linkage group. As an strategy to improve marker ordering, we will find where the major gaps are, and fix it by moving the block of markers to its most likely position when looking at the heatmap.
Linkage group 1
From the map summary tables, we notice that from the first run of Branch-and-Bound (maize.bb1) provided a better map than the second run (-2697>-2711 ) even though its length was slightly greater (215>214 ).
plotRF(maize.bb1, chr=1, col.scheme = "redblue")
By looking at its heatmap, there seems to be small inversions, but no bigger inversions that can be easily corrected manually, so we keep it as is. Small inversions can be fixed by the function ripple(). However, this function is not well optimized in R/qtl for many markers, so we can avoid it at this point.
Linkage group 2
For linkage group 2, the first run also provided greater log-likelihood (-2006>-2016), with a length of 259cM and a gap as big as 91 cM, so we need to fix it. From the heatmap, we notice that there are three subgroups of markers, where the one in the middle seems misplaced. We locate the markers on each side by looking at the major gaps.
plotRF(maize.bb1, chr = 2, col.scheme = "redblue")
(gaps <- tail(sort(diff(maize.bb1$geno[[2]]$map)), 2))## PZA02480.1 PZA03317.1
## 14.60649 18.12182match(names(gaps), names(maize.bb1$geno[[2]]$map))## [1] 137 106maize.bb1 <- switch.order(maize.bb1, chr = 2, order = c(1:10,35:139,34:11), error.prob = 0)pull.map(maize.bb1, chr = 2)## PZA01887.1 PZA00818.1 PZA01983.1 PZA02367.1 PZA00191.5
## 0.00000000 0.00000005 0.00000010 0.00000015 5.19673118
## PZA01438.1 PHM5359.10 PZA02316.22 PZA01570.1 PHM13122.43
## 8.78602208 10.12026021 10.12026026 12.24409378 13.40604097
## PZA00934.2 PZA00801.1 PZA02113.1 PZB00869.4 PZA00499.3
## 66.62887524 68.07610516 68.07610521 68.07610526 68.42002535
## PHM16854.3 PZA01563.1 PHM4647.8 PHM12992.5 PZA00981.3
## 68.42002540 68.42002545 68.42002550 68.73594392 68.73594397
## PZA00996.1 PHM3691.18 PZA01530.1 PZA02207.1 PZA00805.1
## 68.73594402 68.73594407 68.73594412 71.46872414 74.51439510
## PZA00522.12.7 PZA01804.1 PZA00222.7 PHM3171.5 PZB01115.3
## 74.83888015 74.83888020 75.12460994 76.60525407 76.60525412
## PZB00232.2 PZA03451.5 PZA02676.2 PHM1870.20 PZA01050.1
## 76.60525417 76.89098391 76.89098396 76.89098401 77.20032273
## PHM4165.14 PZB01112.1 PZA01349.2 PZA03677.1 PZA02525.1
## 77.56009639 78.66821768 79.02829949 79.38838130 79.38838135
## PZA01303.1 PHM5798.39 PZA00261.6 PZA02862.3 PZA02818.6
## 79.38838140 79.38838145 80.87189536 80.87189541 81.45799619
## PZA01779.1 PZA00273.5 PZA01693.1 PZA03049.24 PZA00643.13
## 81.97014541 83.32644453 83.32644458 84.68274369 84.98422789
## PZA02164.16 PZA00881.1 PZA01365.1 PZA00067.10 PZB01017.1
## 85.59297131 85.59297136 88.20603268 90.90651108 94.51805615
## PZA01608.1 PZA02981.2 PZA01796.1 ae1.8.7 PZA00148.3
## 94.51805620 95.79137298 95.79137303 98.27741625 98.27741630
## PZA01763.2 PZA03714.1 PZA02209.2 PZA02408.2 PZA01575.1
## 98.27741635 104.50469474 105.72315886 105.72315891 106.42507763
## PZA03452.6 PZA02356.7 PHM5296.6 PZA03324.1 PZA02633.4
## 107.86543608 110.13819492 110.13819497 110.13819502 110.50330329
## PHM1899.157 PZA01304.1 PZA03717.1 PZA00987.1 PZA00300.14
## 111.15128634 111.15128639 112.61126785 114.88989046 115.42436092
## PZA02040.2 PZA01410.1 PZA00255.14 PZA01294.2.1 PZA03536.1
## 115.42436097 115.95883142 117.33041552 117.61386222 117.61386227
## PZA02641.2 PZA03317.1 PZA03172.3 PZA02383.1 PZA02751.1
## 117.61386232 138.76026121 139.36437435 139.76017155 140.60423825
## PZA00352.23 PZA03320.6 PZA02426.1 PZA02411.3 PHM532.23
## 140.60423830 140.60423835 142.56091706 143.36480135 144.16868564
## PZA03024.16 PZA00652.17 PZA01142.4 PZA01265.1 PZA02820.17
## 148.94221112 148.94221117 149.27240851 151.95729104 151.95729109
## PZA02513.1 PZA02060.1 PZB00765.1 PZA02667.1 PZA00395.2
## 152.25251498 152.25251503 152.25251508 152.89975815 154.08609088
## PHM3512.186 PZA00963.3 PZA00980.1 PZA02068.1 PZA00545.26
## 154.86526362 156.63975793 156.97466917 158.49701457 165.56584623
## PZA02015.11 PZA03167.5 PZA01259.1 PZA00836.1 PZA01680.3
## 167.60771018 169.08468178 170.17733327 170.50537685 170.50537690
## PZA01140.1 PZA02390.1 PZA02480.1 PZA02769.1 PZA01060.1
## 178.60589687 179.27267561 196.04772402 197.59857709 197.59857714
## PHM565.31 PZA01427.1 PZA02792.26.25 PZA03274.4 PZA03226.3
## 312.90195814 312.90195819 314.77629044 316.65062269 316.94138978
## PZA00517.7 PZA01523.1 PZA03578.1 PZA01327.1 PZA00985.1
## 318.48717117 319.99698330 320.48880292 321.48406631 324.85894197
## PZA03092.7 PZA00112.5 PZA01284.6 PZB00079.4 PZA00865.1
## 326.67015311 326.67015316 331.51131503 333.75956833 335.83115859
## PZA01925.1 PZA01371.1 PZA02029.21 PHM3137.17 PZB00094.1
## 337.90274886 337.90274891 337.90274896 338.20743803 338.20743808
## PZA02753.1 PZB00054.3 PZA02462.1 PZA02653.12
## 340.24575335 341.48102964 343.72285027 345.96467091plotRF(maize.bb1, chr = 2, col.scheme = "redblue")
(loglik.bb1[2] <- attr(maize.bb1$geno[[2]]$map, "loglik"))## [1] -2469.782And we can keep doing this for every linkage group. However, another alternative delivers better ordering in general much more efficiently.
Multidimensional Scale (MDS) by MDSMap + MAPpoly
MDSMap implements the MDS algorithm, and MAPpoly has a nice wrapper function called mds_mappoly. This function only needs the recombination fraction matrix estimated by R/qtl, and it returns a new order that can be used by switch.order() R/qtl function to yield a newly ordered linkage group.
First, we need to install MAPpoly package and use the function below to get the order given by the MDS algorithm for a chosen linkage group. You need to run it first so it becomes available into R before you actually use it:
library(mappoly)## Warning: package 'mappoly' was built under R version 4.1.3## =====================================================
## MAPpoly Package [Version 0.3.02022-01-11 09:22:41 UTC]
## More information: https://github.com/mmollina/MAPpoly
## =====================================================getMDSorder <- function(cross, chr){
markers <- match(names(cross$geno[[chr]]$map), colnames(cross$rf))
mat <- cross$rf[markers,markers]
rec.mat <- lod.mat <- matrix(rep(NA, length(markers)^2), nrow =length(markers))
colnames(rec.mat)<-colnames(lod.mat)<-rownames(rec.mat)<-rownames(lod.mat)<-colnames(mat)
lod.mat[upper.tri(lod.mat)] <- mat[upper.tri(mat)]
lod.mat[lower.tri(lod.mat)] <- t(lod.mat)[lower.tri(lod.mat)]#; image(lod.mat)
rec.mat[lower.tri(rec.mat)] <- mat[lower.tri(mat)]
rec.mat[upper.tri(rec.mat)] <- t(rec.mat)[upper.tri(rec.mat)]#; image(rec.mat)
input.mat <- NULL
input.mat$rec.mat <- rec.mat
input.mat$lod.mat <- lod.mat
mds.map <- mappoly::mds_mappoly(input.mat)
mds.ord <- match(as.character(mds.map$locimap$locus), colnames(mat))
return(mds.ord)
}We’ll create a new object called maize.mds which is a copy of our original cross object maize , so that we can update the ordering within maize.mds only. In addition, we’ll create an empty object called loglik.mds to store the log-likelihood of the orderings obtained using MDS:
maize.mds <- maize
loglik.mds <- c()Linkage groups 1-10
for(i in 1:10) {
c <- i
mds.ord <- getMDSorder(cross = maize.mds, chr = c)
maize.mds <- switch.order(cross = maize.mds, chr=c, order=mds.ord, maxit = 10000, tol=1e-5)
plotRF(maize.mds, chr=c, col.scheme = "redblue")
pull.map(maize.mds, chr = c)
(loglik.mds[c] <- attr(maize.mds$geno[[c]]$map, "loglik"))}## Stress: 0.63686
## Mean Nearest Neighbour Fit: 3.6199
## Stress: 0.14514
## Mean Nearest Neighbour Fit: 3.85791
## Stress: 0.18844
## Mean Nearest Neighbour Fit: 3.42242
## Stress: 0.72087
## Mean Nearest Neighbour Fit: 4.85077
## Stress: 0.63882
## Mean Nearest Neighbour Fit: 5.48484
## Stress: 0.60865
## Mean Nearest Neighbour Fit: 4.29547
## Stress: 0.22367
## Mean Nearest Neighbour Fit: 4.26468
## Stress: 0.27102
## Mean Nearest Neighbour Fit: 4.71302
## Stress: 0.63129
## Mean Nearest Neighbour Fit: 5.63006
## Stress: 0.35705
## Mean Nearest Neighbour Fit: 4.1878
pull.map(maize.mds, chr = c)## PZA02554.1 PZA02221.20 PZA01313.2 PHM3631.47 PHM2828.83 PZA02095.10
## 0.000000 0.266806 7.426917 10.349886 26.111106 26.111106
## PHM3765.7 PZA01883.2 PZA01451.1 PZB01301.5 PHM15331.16 PZA00463.3
## 26.111107 27.642752 27.642753 31.404038 31.404038 37.757377
## PHM3896.9 PZA01642.1 PZA02961.6 PHM3922.32 PZA00079.1 PZA02470.2
## 39.634301 39.634302 39.634302 39.634303 40.679807 40.679808
## PZA02853.11 PZD00033.3 PZA03491.1 PZA00409.17 PZA01597.1 PHM2770.19
## 44.327825 44.690121 45.050062 45.398393 45.398394 46.134476
## PZA01677.1 PZA02941.7 PZA00933.3 PZA01877.2 PZB00409.6 PZA00337.4
## 46.134476 46.134477 46.134477 48.415363 51.961892 52.779302
## PHM12990.15 PZA00814.1 PZA01619.1 PZA00048.1 PZA00400.3 PZA02398.2
## 52.779302 53.130049 53.496857 53.496858 53.496858 55.785580
## PHM537.22 PHM12625.18 PHM18195.6 PZA00444.1 PZA01292.1 PZA01919.2
## 55.785581 55.785581 56.149257 56.149258 56.838156 56.838242
## PZA02128.3 PHM13687.14 PZA02219.2 PZA01089.1 PHM4341.42 PZA01141.1
## 57.145309 61.371671 61.371671 61.371672 61.371672 61.965458
## PZA03196.1 PZA03713.1 PZA01005.1 PZA00866.2 PZA00647.9 PZA01241.2
## 66.093122 66.093122 66.093123 66.093123 66.093124 67.063494
## PZA02320.1 PZB01111.8 PZA01456.2 PHM18513.156 PHM15868.56 PZA02663.1
## 67.063495 71.524491 73.405765 73.962572 74.951572 81.500763
## PZA01995.2 PZA03603.1 PZA03605.1 PZA03606.1 PZA03604.1 PZA03607.1
## 81.500764 86.255385 86.255385 86.255386 86.255386 86.255387
## PZA02049.1 PZA02969.9 PZA00130.9 PZA00007.1 PHM5435.25 PZA01073.1
## 86.255387 86.255388 86.255388 86.255389 95.058912 96.645705
## PZA02167.2 PZA01001.2 PZA02578.1 PZA00062.4 PZA02527.2
## 99.468548 106.541547 106.825264 109.653371 115.607264loglik.mds[c] <- attr(maize.mds$geno[[c]]$map, "loglik")knitr::kable(cbind(summaryMap(maize.mds), log.likelihood=c(loglik.mds, sum(loglik.mds))))| n.mar | length | ave.spacing | max.spacing | log.likelihood | |
|---|---|---|---|---|---|
| 1 | 175 | 1617.9053 | 9.298306 | 171.77421 | -5043.518 |
| 2 | 139 | 166.8328 | 1.208933 | 16.71022 | -2177.542 |
| 3 | 130 | 286.9947 | 2.224765 | 109.02960 | -2343.947 |
| 4 | 127 | 219.5583 | 1.742526 | 31.33065 | -2442.516 |
| 5 | 111 | 162.1694 | 1.474267 | 8.75744 | -2077.314 |
| 6 | 106 | 143.3530 | 1.365266 | 15.47375 | -1694.071 |
| 7 | 85 | 129.0943 | 1.536837 | 20.52314 | -1489.039 |
| 8 | 78 | 123.8148 | 1.607985 | 12.05163 | -1395.635 |
| 9 | 78 | 151.3986 | 1.966215 | 17.29569 | -1517.055 |
| 10 | 77 | 115.6073 | 1.521148 | 15.76122 | -1400.685 |
| overall | 1106 | 3116.7283 | 2.843730 | 171.77421 | -21581.323 |
plotMap(maize.mds)
plotRF(maize.mds, col.scheme = "redblue")